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The present invention is concerned with matched filtering. One example of the use of matched 
filtering is in imaging applications, where an original signal is input to some system which 
produces a response and one wishes to characterise the system by comparing the response with the 
original signal. Some examples of imaging systems are telecommunications line testing and optical 
time domain reflectometry, where the original signal is transmitted to the line, and the response 
obtained from the far (or near) end of the line. Another is ranging systems such as RADAR or 
SONAR where the original signal is transmitted as a radio or ultrasonic wave and the response 
consists of reflections of it from the surroundings. Others include magnetic recording channel 
characterisation, and spatial separation (spectroscopy). 

In general, the response may include the original signal, echoes of the original signal, distortion 
and noise. Conceptually, the characterisation involves examining the correlation between the 
original signal and the response in order to distinguish those aspects of the response that are a 
result of the original signal, rather than the result of some other signal, or other interference or 
noise. Noting that the correlation (or convolution) of two signals in the time domain is equivalent 
to a multiplication operation in the frequency domain, one can perform this task by multiplying the 
frequency spectrum of the response by the frequency spectrum of the original signal A common 
approach to implementation of matched filtering in this way is to form the Fourier transform of the 
response, multiply it by the Fourier transform of the original signal, and take the inverse Fourier 
transform of the product. A particularly attractive approach is to deal with signals containing a 
number of signal samples equal to a power of two, so that the transforms may be performed using 
the fast Fourier transform (FFT). 

Of particular interest is the situation where the original signal consists of a series of pulses. In the 
analysis given below, such a sequence is considered as a series of Dirac impulses (i.e. idealised 
pulses of infinitesimal duration) of magnitude +1 or -1 at regular time intervals. In practice, of 
course the original signal will consist of a band-limited version of such an idealised sequence. For 
convenience of notation, such a sequence will sometimes be expressed as a binary number -where, 
for example, a pulse sequence (1, 1, 1,-1) would be represented as 1 1 10. 

The Golay complementary sequences are well known (see MJ. E. Golay, "Complementary 
Series", IRE Transactions on Information Theory, vol. 7, pp. 82-87, April 1961). They are pairs of 
finite binary sequences with certain useful autocorrelation properties. 
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We define the following notation: 

X 4 " is a binary sequence (e.g. 1110); X" is its inverse (e.g. 0001). 

The concatenation of two such symbols represents the concatenation of the binary sequences - thus 
5 X*X" would, with the example values just given, represent 1 1 100001 . 
x(t) is the time domain sequence corresponding to X* (e.g.(+l,+l,+ 1,-1). 
X(f) is the Fourier transform of x(t). 

A Golay pair consists of two sequences A + and B + , each of length N, having the property that the 
sum of the autocorrelation of a(t) and the autocorrelation of b(t) (both computed at a shift k) is 2N 
1 0 for k=0 and zero for all |k|>0. Examples of Golay pairs having lengths equal to a power of two 
may be generated by the following iteration formula (written in the binary notation), where a pair 
in which both sequences have 2 K bits is written (A + k,B + k ): 

code of length 2° (=1): let (A + 0 ,B + 0 ) = (1,1) 

code of length 2 K+1 : let (A + k+ i,B + k-h) = (A + K B + K , A + k B'k) 

15 For example, (A + 0 ,B + 0 ) = (1,1) 

(A + l9 B + 0 = (11,10) 
(A + 2 ,B + 2 ) = (1110,1101) 
(A + 3 ,B + 3 ) = (11101101, 11100010) 

etc. 

20 This iteration produces only one pair for a particular value of K, whereas, in general, there are 

more: other examples can be generated from these examples using the procedures described in the 
Golay paper cited above. 

According to the present invention there is provided .... 

Some embodiments of the present invention will now be described, by way of example, with 
25 reference to the accompanying drawings. 

Conventionally, matched filtering may be performed in a number of ways, for example by means 
of a finite impulse response filter having tap weights chosen to correspond to the reference signal, 
or using the Fourier transform. Figure 1 is a block diagram of a conventional matched filter of the 
latter type. 

30 A signal r(t) to be filtered (perhaps the response received in an imaging application), and sampled 
at intervals x is received at an input 1 and subjected at 2 to the fast Fourier transform to produce a 
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frequency domain signal R(f). A reference signal s(t) (perhaps the original signal in an imaging 
application) representing the signal that the filter is to match, and consisting of regular pulses at 
intervals X, is received at an input 3 and subjected at 4 to the fast Fourier transform to produce a 
frequency domain signal S(f). Note that the different symbols x, X are used for generality; in the 
5 conventional system X=x and the Fourier transforms are generated for discrete frequencies which 
are multiples of l/2N?i= l/2Nx. It is observed that, in the imaging application, this is perfectly 
satisfactory provided that the receiver of R(t) has a timing reference that is synchronised to that 
used to generate the original signal s(t), and the phenomena being observed are perfectly 
stationary. 

1 0 The product Y(f) = R(f)S(f) is formed by multipliers 5 and the product is transformed by an inverse 
fast Fourier transform at 6 back to a time domain signal y(t) which is output at 7. 

The embodiment of the present invention now to be described is a matched filter of construction 
similar to that of Figure 1, but avoids the use of the FFT in box 4. Moreover it admits of the 
possibility that the sampling period x used for the signal r(t) to be filtered is not the same at the 
1 5 pulse period X of the reference signal s(t). This can be useful in imaging applications because it 
removes the necessity for synchronisation of the analogue parts of the receiver, since small timing 
differences can be accommodated digitally at a later stage in the processing; even with 
synchronisation, it can be useful for accommodating apparent timing differences due, for example 
due to Doppler shifts (as in RADAR). 

20 Thus, in Figure 2, the FFT 4 is replaced by a computation unit 14 which synthesises the Fourier 
transform of the reference signal s(t). Note that in this case the reference signal s(t) does not 
explicitly exist as a time domain signal; the computation only requires its binary representation S + 
(note that the notation S + is used to distinguish it from the Fourier transform S(t): S + is not 
necessarily one of a Golay pair) . 

25 Preparatory to explaining how this computation is performed, we revisit the question of the 

structure of the pulse sequence. The sequence consists of a series of Dirac impulses spaced at time 
intervals of X, with a bit = 1 represented by an impulse of amplitude +1 and a bit = 0 represented 
by an impulse of amplitude -1. An impulse att=0 is written 5(t): whilst an impulse at time T is 
written 8(t-T). (The term "amplitude " is used for convenience: strictly, a Dirac impulse has 

30 infinite amplitude and infinitesimal duration, but the product of these is finite. References to unit 
amplitude refer to a pulse in which this product is equal to unity). 

For a finite bit sequence the impulse train is by default centred on time t= 0, so, for example: 
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binary form Function of time 

This example sequence is shown in Figure 3. 

Any binary codeword can be generate from a single "1" by a combination of concatenation and 
inversion steps (and of course this also applies to the iteration formula given above for generating 
examples of Golay sequence pairs). 

5 Bitstring inversion can be achieved in the equivalent time domain signal by multiplying the time 
function by -1. In the frequency domain it can be achieved by multiplying the Fourier transform 
by -1. This if a binary sequence X* is represented in the time domain by x(t) having Fourier 
transform X(f) then the inverse X" is represented in the time domain by -x(t) having Fourier 
transform -X(f) 

1 0 Concatenation of two bitstrings X* and Y* can be achieved in the time domain by shifting the first 
function back in time, the second function forward in time, and adding them. If the first string has 
N x bits and the second has N Y bits then the shifts are A,N x /2 and A,N Y /2 respectively. In the 
frequency domain one multiplies the first spectrum by exp(+j27if \N x /2) and the second by exp(- 
j27cf XNy/2) and then adds them. 

1 5 Now, the Fourier transform of a single impulse 8(t) is 1 . 



The above equivalents are tabulated in Table 1 : 



operation in 


bitstring domain 


time domain 


frequency domain 


notation 


X^ 


x(t) 


X(f) 


initially 




5(t) 


Kf) 


invertCX*) 


x- 


-x(t) 


-X(f) 


concatenate( A + , B* ) 
where A + has N A bits 
where B + has N B bits 


A + B + 


a(t+^N B /2) 
+ 

b(t-XN A /2) 


A(f).exp(+j7tA.N B f) 
+ 

B(f).exp(-jrt^N A f) 
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Using these one can derive an iteration formula for the Fourier transforms of the Golay sequence 
pairs equivalent to that give earlier for the binary representation: 

setA 0 (f):=1.0+jO.O 

set B 0 (f) := 1.0 +j 0.0 
5 set <|>o := V* X 

(<{> being the period to shift the time functions in the next concatenation) 

forK:=ltoNdo 

set A K (f) := A K -i(f) exp(+j27c<|>K-if) + B K -i(f) exp(-j27t(|>K-if) 
set B K (f) := A K -i(f) exp(+j27c<|>K-if) - B K -i(f) exp(-j27u}>K-if) 
10 set()) K :=2<()K-i 

This computation needs to be performed for each value of f. For a sampled sequence the frequency 
spectrum is theoretically infinite, but repeats. Any frequency interval of width twice the Nyquist 
frequency can be taken, though conventionally it is usual to evaluate f = n/2xM (n = -M+l ...M) 
1 5 where M is the number of samples embraced by the FFT analysis at 2 (Figure 2). Note that it is not 
necessary to evaluate for n = -M since the result is the same as for n =M; and that - noting that the 
transform for -f is the complex conjugate of that for f, one may choose to perform the recursion 
only for positive values of f and then derive those for negative values by taking the complex 
conjugate, 

20 Note that there are alternative versions of the Fourier transform which are scaled versions of one 
another. The usual convention for FFT processing is that the Fourier transform of the unit impulse 
is 1/M; if this is desired then A 0 (f) and B 0 (f) should be initialised to 1/M rather than unity. Note at 
this point that we prefer, in the test method discussed, so send an effectively infinite sequence 
consisting of the concatenation A + B + A + B" sent repeatedly (by "effectively infinite" we mean that 

25 transmission commences a sufficiently long time before measurement begins that transient effects 
have died away, and does not cease until after measurement is complete). (We prefer this 
sequence as it has the property of having no sidelobes within 1/4 cycle of t=0). Lest it should 
appear surprising that the sequence is expressed in this form, a word of explanation is in order. We 
observe that, if the recursion formula given above is applied to the pair (A P + ,B P + ), of some length 

30 2 P , to produce a pair of twice the length (Ap+i^Bp+i*), then transmission of the concatenation 
Ap+i^p+i* amounts to the same thing as sending Ap^p^A.p'lBp-. Moreover, if one applies the 
recursion formula again to (A P+ i + ,B P+ i + ) to give a pair (A p+2 + 3p+2 + ) of twice the length, then 



WO 2004/088530 



PCT/GB2004/000954 



transmission of Ap+ 2 + alone also amounts to the same thing as sending A P + B P + Ap 4 B P ". However, 
whilst all pairs generated by this particular recursion have the property that if one takes one 
sequence of the pair and cuts it in half then one has a new (shorter) Golay pair, it is not in general 
true that all Golay pairs of length 2 n have this property (or at least, it has yet to be proved that this 
5 is so). Thus it cannot be assumed that a single Golay sequence will necessarily have this property. 
However this property can be guaranteed providing that the recursion is performed at least twice. 
Thus if an arbitrary Golay pair of any length, which may or may not have this property, is 
subjected twice to the recursion, the result will have this property. Except where the starting pair 
are of unit length, it follows that (given that the length of a member of a Golay sequence pair is 
1 0 always one or an even number) each member of the resultant pair will have a length of 8 or a 
multiple thereof. In this case the transmitted sequence is Ap+2 + A P+2 + Ap+2 + etc. The sequence to 
which the filter is to be matched is the sequence Af^B^Af^Bf": i.e. in this instance simply A P+2 + . 
So the computation at 14 is required only to compute the Fourier transform Ap+ 2 (f). 

As note earlier, signals occurring in reality do not consist of ideal Dirac impulses; any transmitted 
15. signal, nominally a sequence of such pulses is in fact a sequence of pulses each having a finite 
amplitude and frequency spectrum. Such a signal can be viewed mathematically as the 
convolution of the idealised Dirac sequence convolved with the shape pj(t) of a single pulse - and, 
if need be also with a pulse shape p R (t) to take account of pulse-shaping that may occur at the 
receiving end. In the frequency domain, the convolution of a(t) and pr(t) becomes simply the 
20 product of A(f) and P(f), where Pr(f) is the Fourier transform of pr(t). It is standard practice to 
take this into account in any matched filtering process, but, for the sake of completeness, it is 
pointed out that, in Figure 2, this may readily be implemented by multiplying S(f) by the complex 
conjugate of the Fourier transform of the pulse shape before applying it to the multiplier 5, i.e. 

S(f)Px(f)P R (f). 

25 It has already been noted that the recursion process produces only one example of a Golay pair of 
any particular length. For many applications this is sufficient - i.e. if one sequence of length 2 K is 
as good as any other then there is no motivation to seek others. However it may be that an 
alternative sequence is desired. For example one may want to transmit two sequences over the 
same system and be able to distinguish between them upon receipt. In such a case (as well as the 

30 option mentioned earlier of choosing different starting points for the recursion, Golay in his above 
mentioned paper offers six (time-domain) transformations which may be used singly or in 
combination to generated alternative sequences. 
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These too may be implemented in the frequency domain. Supposing that one has generated a 
Golay pair of Fourier transforms A(f),B(f) equivalent to a time domain pair a(f),b(f), then one can 
generate new pairs A'(f)3'(f) 5 equivalent to new time-domain pairs a»(t),b'(t) as indicated in the 
following table: 





time domain as proposed by Golay 


frequency domain 


a) 


interchange: a'(t) = b(t); b'(t) = a(t) 


A'(f) = B(f);B'(i) = A(f) 


b) 


time reversal of a(t) 


A(f) = A*(f) (where * indicates the complex 
conjugate) 


c) 


time reversal of b(t) 


B'(f) = B*(f) n 


d) 


inversion: a'(t) = -a(t) 


A'(f) = -A(f) 


e) 


inversion: b'(t) = -b(t) 


B'(f) = -B(f) 


f) 


invert alternate bits of both a and b, that is 
to say a'(t) - a(f).g(t) and b'(t) = b(t).g(t) 
where g(t) = (-1,1,-1,1,-1,1 etc.) 


convolve A(f) with G(f) (& similarly for B). 
Since G(t) is a sinusoid of frequency 1/2X, 
its Fourier transform G(f) is simply two 
Dirac impulses at ± f. Thus A(f)*G(f) is 
given by ± JAtf - T JA(f + . 



5 



Some variations on this process will now be discussed: 

The starting point for the recursion above was (Ao(f),Bo(f)) = (1,1). In fact one can start with any 
Golay pair, of any length. Inter alia this includes (Ao(f),B 0 (f)) = (1,1) or (1,-1) or (-1,1) (-1,-1). 

Windowing in t: Usually we shall be interested in analysing data, obtained by measurement, to 
identify instances of a codeword waveform somewhere in it. A naive approach would be to use a 
code cycle which just fits the FFT; this may even be a good idea for applications where the source 
and analyser can be sure of identical clock frequencies. However for the general case we will want 
a guard period at the beginning and end of the time block to prevent spurious responses from 
partial codewords. This may be achieved by using a block which is a multiple (at least twice) of 
the length of the codeword being sought, and analysing with a matched filter representing a single 
cycle and zero outside that cycle (the analyser cycle is assumed centred at time t = 0). The guard is 
achieved by ignoring any response within half a cycle of the block edges. 



10 



15 
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Windowing in f: The shift operations, and their exp(-j2rc shift f) multiplications, will in general 
produce discontinuities in the spectrum at the Nyquist band edges. This will not happen when the 
X value is an exact multiple of t (when the exp( ) factors join smoothly at the edges) but we will in 
general be interested in cases where X is only approximately equal to x. 

5 As we are interested here in real valued signals in time, we would expect Hermite symmetry to be 
maintained in any spectrum. In particular we would expect X[0] and X[f Nyquist ] to be real. So at the 
very least we would expect the values stored at the Nyquist frequency to be forced to be real, 
whatever the process actually calculates. 

The usual solution for such a discontinuity is windowing in the frequency domain - i.e. tapering 
10 the spectrum to zero as the band edges are approached. The consequence of not doing so, and just 
tolerating the discontinuity at f Ny qui s t is that any narrow time domain feature will appear as a sinc() 
function displaced from the time sampling grid - it will acquire a skirt of sidelobes, providing a 
slowly decaying noise floor at about -15dB below the feature's peak. 

There are many suitable window functions in the literature, all with slightly different virtues. The 
1 5 root raised cosine window is a convenient one to choose. 

Golay' s paper also discusses alternative recursions for the generation of Golay sequences. For 
example, the recursion may, rather than concatenating the bits of A + and B + (or B*), interleave the 
bits. In the time domain, this amounts to expanding the pulse sequence along the time axis, 
shifting, and adding. This scheme too may be implemented by analogous operations in the 

20 frequency domain, though it may be simplified by, rather than regarding the time domain operation 
as taking a sequence (for example) A + and B + with constant pulse spacing and successively 
doubling the sequence period, one can instead regard it as generating sequences of constant period 
and pulse spacing halving at each step. If starting with a sequence of length 1 and wanting a Golay 
pair each of 2 N pulses, so there are to be N iterations, then one desires a sequence period of 2 N X 9 

25 and starts with a notional of pulse spacing of 2™X; the first shift is 2 N A/4. 

The recursion is as follows: 

setA 0 (f):= 1.0 +j 0.0 
set B 0 (f) := 1.0 +j 0.0 
set 0 O := 2 N " 2 X 

30 (0 being the period to shift the time functions in the next concatenation) 



WO 2004/088530 



PCT/GB2004/000954 



" 9 " 

forK:= 1 toNdo 

set A K (f) := AK-i(f) exp(+j27te K -if) + B K -i(f) exp(-j2rcG K -if) 
set B K (f) := A K _i(f) exp(+j27te K -if) - B K -i(f) expH27c9 K -if) 
set 0 K := 6 K _,/2 

5 

As before, one may initially set Ao(f),Bo(0 to be the Fourier transforms of any Golay pair. Note 
however that this process will not necessarily produce sequences having the sidelobe property 
discussed earlier, and if this is required this can be achieved by generating sequences which are a 
quarter of the desired length and subsequently performing two iterations of the recursion described 
1 0 earlier. 
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CLAIMS 

1. A method of matched filtering comprising: 
computing the Fourier transform of a signal to be filtered; 

5 computing the Fourier transform of a reference sequence to which the filter is to be matched; and 
forming the product of the two transforms; 

characterised in that the reference sequence is defined by a process of iteratively combining shifted 
versions of shorter sequences and the step of computing the Fourier transform of the reference 
sequence comprises an iterative process of combining the Fourier transforms of a shorter starting 
1 0 sequence. 

2. A method according to claim 1 in which the reference sequence is a Golay sequence pair and the 
step of forming the Fourier transform of the reference sequence comprises repeatedly: 

(a) combining the Fourier transform of a first member of a Golay pair with the Fourier transform of 
1 5 the second member of that Golay pair to produce a first member of a new Golay pair; and 

(b) combining the Fourier transform of a first member of a Golay pair with the Fourier transform 
of the second member of that Golay pair to produce a second member of a new Golay pair. 

3. A method according to claim 2 in which said combining uses only the operations of inverting, 
20 addition, and multiplication by exp(±j27if<j>), where f is frequency and § is a shift value dependent 

on the length of the sequence. 

4. A method according to claim 3 in which the transforms A K (f), B K (f) of a Golay pair are formed - 
form the transforms A K -i(f), B K -i(f) of a shorter such pair according to the relationships 

25 A K (f) := A K _!(f) exp(+j27i(|>f) + B K -i(f) exp(-j27i<|>f) 
B K (f) :« A K -i(f) exp(+j27r([)f) - B K _i(f) exp(-j27i<j)f) 

where $ is half the length of each member of the shorter pair, and f is frequency . 

5. A method according to claim 3 in which the transforms A K (f), B K (f) of a Golay pair are formed 
30 from the transforms A K _,(f), B K -i(f) of a shorter such pair according to the relationships 

A K (f) := A K -i(f) exp(+j27c0 K -if) + B K -i(f) exp(-j27i;0K-,f) 
B K (f) := A K -i(f) exp(+j27c9 K -if) - B K -i(f) exp(-j27c9 K -if) 

where the 9 are time intervals dependent on the number of iterations, and f is frequency . 
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6. A method according to claim 4 or 5 hi which the iteration commences with a Golay pair each 
member of which has a length of 1 . 

5 7. A method according to claim 2 in which said combining uses only the operations of inverting, 
addition, and multiplication by exp(±j2ftf$), where f is frequency and $ is a shift value dependent 
on the length of the sequence. 

S. A method according to claim 7 in which the transforms A*(f)» B K (f) are formed form the 
1 0 transforms A K -i(0> B K -iO of a shorter such pair according to the relationships 
Ak(0 A K -i(f> exp(+j2*<frl) + B K -i(0 exp0-j2raj>f> 
Bit(f) v Aic-i(f) exp(+j2iw}>0 - B k -i(0 exp(-j2jwtf) 

where <j> is half the length of each member of the shorter pair, and f is frequency 

15 9. A method according to claim 7 in which the transforms A*(fX B K (f) are formed from the 
transforms A K -,(f>* Bie-,(f) of a shorter such pair according to the relationships 
A K (f) :» A K -i(0 exp(-»J2iie K -iO + B*-,(f) exp(-j2nB It - l f) 
BK(f) :- Ak,i(0 exp(+j2itOK-if) - B*_,(f) exp^nS^f) 

where the 8 are time intervals dependent on the number of iterations, and f is frequency . 

20 

1.0. A method according to any one of the preceding claims including the step of forming the 
inverse Fourier transform of the product. 
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